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Abstract 

Model  Predictive  Control  (MPC)  is  the  class  of  control  methods  that  optimize 
a  specified  performance  index  in  order  to  minimize  weighted  future  output  devia¬ 
tions  from  a  setpoint  trajectory.  MPC  operates  on  a  receding  horizon,  calculating 
a  set  of  inputs  at  each  time  step.  The  controller  then  implements  the  first  input 
and  the  process  begins  again.  The  performance  index  can  also  include  a  weighted 
and/or  constrained  control  sequence  which  can  be  of  different  length  than  the  output 
horizon. 

This  thesis  applies  MPC  in  the  inverse  sense  -  known  aircraft  outputs  are  ap¬ 
plied  in  the  performance  index  as  setpoints  in  an  attempt  to  determine  what  control 
histories  caused  those  outputs.  Using  this  method,  aircraft  mishap  investigators 
could  then  have  a  means  of  determining  what  the  control  surface  deflections  were 
throughout  an  incident.  To  accomplish  these  objectives,  MATLAB  and  MATLAB 
routines  axe  used  in  conjunction  with  SIMULINK  to  develop  the  controller  and  sim¬ 
ulate  the  aircraft’s  response.  The  actual  Flight  Data  Recorder  data  from  aircraft 
mishaps  are  utilized  as  proof  of  concept. 
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Solving  an  Inverse  Control  Problem 
using  Predictive  Methods 


/.  Introduction 

1.1  Inverse  Control  Problem 

Recently,  much  interest  has  been  plaeed  in  the  area  of  reconstructing  aircraft 
accidents  based  on  the  information  available  in  the  Cockpit  Voice  Recorder  (CVR) 
and  the  Flight  Data  Recorder  (FDR)  post-crash.  The  CVR  can  provide  audible 
clues  of  the  accident’s  cause  such  as  stall  warnings,  engine  noise,  landing  gear  exten¬ 
sion/retraction,  etc.;  whereas,  the  FDR  provides  the  investigator  with  the  operating 
conditions  leading  up  to  and  during  the  crash  such  as  altitude,  heading,  airspeed, 
etc..  Neither  the  CVR  nor  the  FDR  provide  definitive  information  regarding  the 
control  surface  locations  throughout  the  aircraft  incident.  Aircraft  investigators  cur¬ 
rently  use  the  information  included  in  the  CVR  and  FDR  in  an  attempt  to  determine 
the  cause  of  the  accident  -  be  it  a  control  surface  failure  or  other  cause.  However, 
of  great  value  to  the  investigators  would  be  the  control  histories  which  caused  the 
aircraft  response. 

The  effort  herein  focuses  on  the  particular  application  of  post-processing  air¬ 
craft  data  from  an  incident.  However,  the  derivation  and  algorithm  axe  generic 
enough  so  any  system  could  be  used.  In  fact,  a  distillation  process  on  which  model 
predictive  control  had  previously  been  performed  was  used  in  this  thesis  to  validate 
the  solution  algorithm.  Any  plant  model  of  the  form  discussed  here  could  be  used, 
with  the  methods  described,  to  obtain  the  inputs  from  known  outputs. 
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1.2  Model  Predictive  Control 

This  research  effort  applies  predictive  methods  to  the  inverse  control  problem 
in  an  attempt  to  determine,  analytically,  the  control  surface  deflections  that  caused 
the  aircraft  response  obtained  in  an  FDR.  Model  Predictive  Control  (MPC)  is  such 
a  method,  utilizing  on-line  optimization  of  a  specified  performance  index  in  order 
to  minimize  weighted  future  output  deviations  from  a  setpoint  trajectory  along  a 
prediction  horizon.  Typically,  the  performance  index  also  includes  a  weighted  (and 
possibly  constrained)  control  sequence  across  a  control  horizon  which  can  be,  but 
does  not  have  to  be,  the  same  length  as  the  prediction  horizon.  MPC  has  proven  to 
be  beneficial  as  a  flight  controller  (14),  especially  in  the  realm  of  aircraft  control  in 
the  presence  of  actuator  failures  (2).  This  thesis  details  the  usefulness  of  MPC  for 
the  inverse  control  problem  by  developing  a  performance  index  utilizing  the  FDR 
data  from  an  aircraft  accident  as  the  setpoint  trajectory. 

1.3  Research  Objectives 

The  primary  objective  of  this  thesis  is  to  determine  whether  MPC  is  useful 
in  determining  the  control  histories  that  were  involved  in  an  aircraft  incident.  By 
applying  the  theory  to  be  developed  to  actual  aircraft  crashes  of  unknown  (or  sus¬ 
pected)  cause,  this  objective  will  be  accomplished.  Also  to  be  explored  are  the  effects 
of  including  a  rate  term  in  the  performance  index  for  improved  setpoint  following. 
These  objectives  will  be  accomplished  using  MATLAB  (7)  to  develop  the  controller 
and  SIMULINK  (8)  to  implement  the  controller  on-line  with  the  aircraft  model. 

1.4  Thesis  Overview 

Chapter  2  reviews  the  available  literature  used  in  development  of  this  thesis. 

Chapter  3  details  the  methodology  and  necessary  background  information  for 
the  development  of  an  MPC  controller,  including  modifications  for  system  constraints 
and  control  increments. 
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Chapter  4  delves  into  the  first  application,  a  distillation  column,  for  validation 
purposes. 

Chapter  5  presents  the  applications  dealing  with  two  aircraft  mishaps,  dis¬ 
cussing  the  scenarios,  models  used,  and  simulation  results. 

Chapter  6  offers  conclusions  that  can  be  made  based  on  the  results  of  the  work 
presented  and  suggests  areas  for  future  research  in  the  area  of  inverse  control  using 
predictive  methods. 
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II.  Literature  Review:  Inverse  Simulation 

In  the  context  of  this  thesis,  inverse  simulation  is  the  practice  of  determining  a 
system’s  input  when  only  the  outputs  and  the  system  model  are  known.  Typically, 
inverse  control  is  used  when  the  effort  goes  one  step  further  -  including  feedback 
control  for  disturbance  rejection  and  plant  uncertainty  robustness.  Several  methods 
have  been  studied  in  an  attempt  to  perform  inverse  simulation.  The  most  relevant 
for  the  discussion  here  are  what  have  been  called  the  differential  and  integration 
inverse  methods  and  those  utilizing  optimization. 

2. 1  Differential  Inverse  Method 

The  differential  inverse  method  entails  modeling  the  desired  plant  trajectory  as 
a  set  of  dynamic  constraints  imposed  on  the  equations  of  motion  (1).  This  trajectory 
can  be  defined  analytically,  as  a  given  equation,  or  numerically,  as  a  series  of  specific 
points  which  are  then  smoothed  to  provide  a  continuous  trajectory.  The  constrained 
equations  of  motion  axe  then  integrated  until  the  control  terms  arise  which  can  be 
solved  for  directly. 

The  shortcomings  of  this  method  of  inverse  simulation  deal  with  the  numbers 
of  states,  inputs,  and  outputs  (3).  If  a  typical  plant  is  considered  of  the  form 

x(t)  =  Ax(t)  +  Bu{t)  (1) 

y{t)  =  Cx{t) 

the  solution  for  the  control  to  obtain  the  desired  trajectory,  ud,  resulting  from  this 
method  is 

u{t)  =  B~^[x{t)  —  Aa;(t)]  (2) 

where 

x(0  =  (3) 
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initiates  the  state  derivatives.  The  resulting  input  is  used  with  Equation  (1)  to 
obtain  the  state  derivatives  at  subsequent  time  steps.  Letting  rix  be  the  number  of 
states,  Uy  the  number  of  outputs,  and  n„  the  number  of  inputs,  if  Ux  =  ny  =  n^,  no 
problems  arise  in  the  solution.  However,  if  Ux  >  n„  and  n^.  >  ny,  as  is  typical  of 
aircraft  models.  Equation  (3)  can  no  longer  provide  all  the  state  derivatives  necessary 
to  solve  the  problem.  To  avoid  this  problem,  initial  guesses  of  the  states  are  typically 
made. 

The  differential  inverse  method  has  been  shown  to  provide  solutions  to  the 
problem  at  hand.  However,  it  is  desirable  to  develop  cm  algorithm  that  does  not  con¬ 
tain  shortcomings  of  this  nature  where  the  solution  might  depend  on  initial  guesses 
of  states. 

2.2  Integration  Inverse  Method 

The  integration  inverse  method  develops  an  error  vector  defined  as  the  differ¬ 
ence  between  the  actual  system  output  and  the  desired  system  output  and  attempts 
to  drive  that  error  to  zero  (3).  If 

y{kT)  =  G[u(fcr)]  (4) 

where  G  maps  the  input,  «,  to  the  output,  y,  then  with  the  introduction 
vector,  Fe, 

FeWnikT)]  =  G[u4kT)  -  YoikT)] 
where  Yq  is  the  desired  trajectory.  The  solution  then  becomes 

Un+i{kT)  =  Un(fcr)  -  (J{GK(A:T)]})-' .  FEW^kT)]  (6) 


of  an  error 

(5) 
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where  J{}  is  the  Jacobian  matrix  found  through  partial  derivatives  of  the  output 
approximated  by 

Syi{kT)/Suj{kT)  «  [yi{uj  +  Auj)  -yi(uj)  UtJ/Auj  (7) 

Iteration  on  the  desired  input,  t/n+i,  is  then  performed  until  the  actual  output  and 
the  desired  output  differ  only  by  the  tolerance  required. 

In  this  manner,  the  integration  method  avoids  the  shortcomings  of  the  differen¬ 
tial  method  in  that  the  only  requirement  imposed  upon  the  number  of  states,  inputs, 
and  outputs  is  that  Uy  <  Uu.  Therefore,  a  solution  of  this  type  is  of  more  interest 
for  this  effort. 

Optimization  Methods 

The  final  methods  of  relevance  in  available  literature  are  those  involving  op¬ 
timization.  Similar  to  the  differential  method,  optimization  methods  formulate  the 
solution  by  setting  equality  constraints  on  functions  of  the  state  variables.  However, 
the  solution  is  not  found  by  differentiating  the  equations  but  as  a  general  optimiza¬ 
tion  problem. 

In  this  type  of  solution,  a  cost  function  is  defined  and  is  augmented  by  any  or  all 
of:  constraints  on  initial  and  final  conditions,  path  constraints,  a  dynamic  constraint 
equation,  and  input  control  constraints  (5).  The  necessary  and  sufficient  conditions 
for  optimality  are  then  found  as  well  as  the  input  solution  to  optimally  drive  the 
system  to  the  desired  trajectory.  This  method  does  not  require  time  differentiation  or 
output  derivatives,  and  therefore  the  numerical  difficulties  of  the  previous  methods 
can  be  avoided.  Also  avoided  are  the  sensitivities  of  the  results  to  initial  guesses. 

In  this  thesis,  optimization  is  utilized  via  Model  Predictive  Control,  which 
features  on-line  optimization  of  the  control  inputs  and  real-time  simulation  of  the 
plant  outputs. 
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III.  MPC-State  Space  Formulation 

A  method  of  control  which  incorporates  the  cidvcuitages  of  optimization  and 
eliminates  the  disadvantages  of  time  differentiation  and  initial  guess  sensitivities  is 
Model  Predictive  Control  (MPC).  This  section  presents  the  development  of  a  state- 
space,  constrained  MPC  controller  which  utilizes  control  increments. 

It  should  also  be  noted  that  in  several  instances,  users  of  MPC  have  included 
in  the  problem  a  stabilizing  inner  feedback  loop  (2).  This  is  done  in  order  to  ensure 
that  the  controller  in  these  forward  loop  applications  produces  only  stabilizing  inputs. 
This  will  not  be  the  practice  here  as  application  of  a  stabilizing  inner  loop  might 
prevent  the  algorithm  from  following  an  output  trajectory  which  is  in  fact  unstable. 
Whether  the  setpoint  trajectory  is  stable  or  unstable,  as  could  be  the  case  in  aircraft 
incidents,  the  inverse  controller  should  be  able  to  follow  the  given  outputs. 

In  this  thesis,  the  controller  is  based  on  a  linear  discrete-time  plant  model  of 
the  form 


x{k  +  1)  =  Ax{k) Bu{k)  (8) 

y(^k)  =  Cx{k)  Du{k) 

where  ®  €  .R”,  is  the  state  vector,  u  E  is  the  input  vector,  and  y  £  BP,  is  the 
output  vector  so  A  G  BP^^,B  G  BP^^,  C  G  and  D  G  BP^^. 

The  Model  Predictive  Controller  uses  this  plant  model  (Equation  (8))  to  cal¬ 
culate  the  future  plant  outputs,  y,  along  a  specified  prediction  horizon,  p,  due  to  the 
inputs,  u,  implemented  over  the  control  horizon,  q,  and  minimize  the  cost  function 
involving  the  difference  between  these  outputs  and  the  desired  setpoint  output  tra¬ 
jectory,  yd.  The  outputs  are  weighted  by  Q,  and  the  inputs  are  weighted  by  B  across 
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their  respective  horizons.  This  gives  the  overall  cost  function 

J  =  -£\\Ci(,k  +  l)-y,(k  + 1)  ||«  +  E  II  + 1)  Hr  (9) 

<=1  fco 

Noting  that  the  actual  plant  states,  x,  are  not  assumed  to  be  known  at  any 
time  in  the  future,  it  is  necessary  to  estimate  the  future  states,  x{k  +  €),  based  on 
past  states  and  inputs  and  future  inputs.  Based  on  an  estimator  of  the  form 

x{k  +  1)  =  Axik)  +  Bu{k)  +  L[y{k)  -  Cx{k)]  (10) 


it  is  seen  that 

x{k  +  £)  =  {A-  LC)x{k  +  £-!)  +  Bu{k  +  ^  -  1)  +  Ly{k  +  £-1)  (11) 

With  this  equation,  it  is  then  possible  to  express  the  state  estimates  as  (using  the 
notation  of  (6)) 


x{k-\-^A\)  =  F{£)x{kAt}  +  Gu{k  +  £)  +  H{£)y{k^-£)  (12) 

where  F,  G,  and  H  are  defined  as  F{m)  =  (A  —  XC')’"''"^,  G  =  B^  H(0)  =  —L  and 
H{m)  —  0  for  m  =  !•••£.  These  terms  dealing  with  output  feedback  in  H  are 
eliminated  for  all  times  greater  than  k  +  1  since  the  actual  system  output  will  not 
be  available  at  any  future  time. 

It  is  now  possible  to  define  a  vector  of  future  predicted  states: 


x{k  + 1) 
x{k  +  p) 


=  Fx{k)  +  g[u{kf  ...u{k-\-q-  1)^]^  +  n{£)yik) 


(13) 
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where  JF,  and  H  are  matrix  functions  of  F{€)^  G,  and  H{i)  as  follows: 

n 

(A-L-C) 

{A-L-  cy 

(A-L-  cy 

B  0  0  0 

(A-L -0)8  B  0  0 

{A-L-CfB  {A-L -0)8  0  0  ^pn  (15) 

i  i  0 

{A-L-Cy-^B  {A-L-Cy-^B  •••  B 

_ V 

L 

(A-L’  C)L 
{A-L-  cy-^L 

The  state  prediction  vector  can  then  be  used  to  transform  the  performance 
indices  developed  in  the  next  sections  into  the  form  of  a  quadratic  program 

J  =7  ^{U^PU  +  fU]  (17) 

subject  to  AU  <  b 

where 

U  =  +  1)^  •  •  •  u(A:  +  q  -  1)^]^  (18) 

since  the  index  is  minimized  over  the  control  input,  U.  Also,  P  and  /  will  be  functions 
of  the  state  prediction  matrices.  In  this  form,  MPC  is  easily  implementable  using  a 
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quadratic  program  algorithm  to  determine  the  future  inputs  to  drive  the  system  to 
the  desired  trajectory. 

3.1  Non-rate  performance  index 

Using  the  state  prediction  matrices,  Equation  (13),  and  expanding  the  cost 
function.  Equation  (9),  into  the  form  of  Equation  (17)  results  in 

J  =  u(kf{g'^c^Qcg  +  ii)U(k)  +  2{[«(A:)  +  ny(k)]^c^  -  s}Qcgu{k)  +  «  (i9) 

and  K  includes  all  terms  from  the  expansion  independent  of  t/,  and  is  therefore 
neglected.  In  addition,  the  matrices 


diag{C  •  •  •  C)  (20) 

diag{Q  •••Q) 
diag{R  •  •  •  R) 
yd{k) 
yd(k  4- 1) 

Vdik  +  p) 


3.2  Rate  performance  index 

It  has  been  suggested  that  inclusion  of  an  error  rate  term,  in  addition  to  the 
previously  discussed  output  error  term,  in  the  performance  index  may  improve  the 
controller’s  tracking  of  the  setpoints.  In  order  to  include  an  error  rate  term,  define 
the  error,  e,  as 

e{k-i-e)  =  Cx(k-\-£) -yd{k-\-e)  (21) 


C  = 

Q  = 
?^  = 

5  = 

are  defined. 
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Then,  define  a  performance  index  similar  to  Equation  (9)  «is 

J  =  x;  II  ^e(k  +  <)  +  Ee(k  + 1)  lla  +  x;  II  + 1)  lli  (22) 

/=!  <=0 

where  the  notation  is  identical  to  that  in  Equation  (9)  with  the  addition  of  Ae,  the 
error  rate  term,  and  E,  which  is  a  weight  placed  on  the  output  error.  By  weighting 
the  output  error  in  this  manner,  the  actual  output  approaches  the  setpoint  arbitrarily 
fast  by  choice  of  E  through  Ae  =  —Ee,  since  the  performance  index  attempts  to 
drive  Ae{k  +  ^)  +  Ee{k  +  £)  to  zero. 

The  state  prediction  equation  will  require  some  manipulation  with  the  inclusion 
of  the  error  rate  term,  Ae(k)  =  e(k)  —  e{k  —  1)  or 

Ae{k  +  £)  =  Cx{k  +  £)  —  yd{k  +  £)  —  Cx{k  +  £  —  !)  +  yd{k  +  ^  —  1)  (23) 

which  can  be  expanded  into  a  vector,  Ae(k),  as 

Ae(fc)=  [Ae(A:+l)3’...Ae(jb  +  p)^]^  (24) 

Ae(A:)  =  CJ^x{k)  +  CGU{k)  +  Cny{k)  -y-  C^x{k)  -  CGU{k)  -  CHy(k)  +  y  (25) 

Using  the  notation  of  Equation  (13)  and  noting  T  is  identical  to  time-shifted 
with  the  first  row  as  zeros  due  to  the  time  shift  between  e(ifc)  and  e(A:  —  1): 

n 

- - , 

0 

{A- L-C) 

f=  {A-L-Cf  (26) 

{A- L-  cy~^ 
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And  G  and  H  are  defined  similarly.  Finally,  we  note 

yd{k  + 1)^ 
yd{k  +  2)^ 

yd{k+pY 

yd{kf 

y=  \ 

yd{k-\-p-\) 

We  can  then  rewrite  Equation  (25)  as 

Ae(A:)  =  C[Tx{k)  +  GU{k)  +  Hyik)]  +  y  (29) 

AAA  A  ^ 

if  we  let  .F ,G,'H,  and  y  equal  T etc.  Also,  from  direct  substitution  of  Equation 
(13)  into  (21),  across  the  horizon,  we  have 

e{k)  =  C[Tx{k)  +  GU{k)  +  Hy{k)]  -  y  (30) 

Utilizing  Equations  (29)  and  (30)  and  expanding  the  cost  function.  Equation 
(22),  into  the  general  quadratic  program  form,  we  see 

J  =  u{kf[g'^c'^Q{CG  +  2ecG)  +  g'^c^s^qscg  +  ii]u{k)  (31) 
-\-2{[CTx{k)  +  CHy(k)  +  yYQ[CG  +  eCG] 

+  [eTx{k)  +  CHy{k)  -  y]^  S^Q[CG  +  SCG]}U{k)  +  k 

and  K  once  again  includes  the  neglected  terms  independent  in  U.  This  formulation 
allows  the  MATLAB  functions  which  will  perform  the  quadratic  program  to  calculate 
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the  original  matrices,  etc.  and  manipulate  those  into  the  required  matrices, 
etc.  for  the  index  including  the  error  rate  term. 


3.3  System  Constraints 

Typically,  the  inputs  determined  from  an  MPC  controller  for  an  aerospace 
system  will  be  constrained  by  rate  and  deflection  limits  based  on  either  servo  or 
physical  limitations  in  both  the  maximum  and  minimum  positions.  This  is  in  or¬ 
der  to  prevent  the  controller  from  attempting  to  drive  the  system  with  inputs  that 
cannot  be  physically  attained.  With  emphasis  on  the  inverse  control  framework,  the 
desire  is  to  reconstruct  inputs  that  have  already  been  imposed  on  an  actual  sys¬ 
tem,  and  therefore  the  inputs  will  only  be  constrained  to  their  physical  deflection 
limits,  both  positive  and  negative.  This  restriction  will  ease  calculations,  decrease 
calculation  time  slightly,  cind  provide  no  adverse  effect  on  the  results  based  on  the 
above  discussion.  In  fact,  cases  can  be  envisioned  where  the  control  deflection  that 
actually  occurred  exceeded  servo  rate  limits  (implying  some  failure),  and  exclusion 
of  those  cases  could  prevent  good  setpoint  tracking.  Also,  the  output  will  only  be 
constrained  at  the  current  time  step  instead  of  across  the  control  horizon  to  allow  the 
controller  to  be  more  aggressive,  while  still  operating  within  limitations,  in  following 
the  potentially  volatile  output  trajectory. 

Finally,  the  constraints  pertinent  for  an  aerospace  system,  written  across  the 
control  horizon,  are  most  easily  expressed  as 


< 

u{k) 

< 

+  5  1) 

u{k  +  q  —  1) 

(J  1) 

(32) 
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These  constraints  must  be  expressed  in  the  same  framework  as  the  quadratic  program 
given  in  Equation  (17)  as  AU  <  6,  or 


I 

-I 


u{k  +  q-  1) 


^maxi^k  ~\~  q  1) 


~Umin{k  +  q  -  1) 


(33) 


Constraints  will  only  be  applied  in  this  thesis  when  working  with  the  trcins- 
port  aircraft  model.  No  constraints  are  placed  on  the  distillation  process  since  it  is 
an  unconstrained  benchmark  problem.  The  physical  deflections  limitations  for  the 
transport,  then,  are  (4) 


•  Elevators  (^e)  =  4-19.33  deg, -21.83  deg 

•  Ailerons  (5o)  =  ±20  deg 

•  Rudder  (6^)  =  ±26  deg 

•  Throttle  (^y)  =  0-4  100% 

3.4  Plant  Modifications  for  Control  Increments 

The  computer  program  used  to  solve  the  quadratic  program  utilizes  an  incre¬ 
mental  control  input  instead  of  absolute  inputs,  so  the  aircraft  model  will  need  to 
be  modified  to  accept  these  inputs.  One  method  of  achieving  this  is  to  introduce  an 
additional  state  corresponding  to  each  input  so 


Xj,{k  ±  1) 

Xj,{k) 

A  A 

Xp(k) 

Xu{k  ±  1) 

=  A 

^ti(^) 

±  BAu{k)  y(k)  =  C 

Xu{k) 
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where 


A  B 
0  I 


B  = 


B 

I 


C  =  [C  0] 


(35) 


In  this  manner,  the  additional  input  states,  x^ik)  =  u{k  —  1),  allow  absolute  in¬ 
puts  to  be  carried  from  one  time  step  to  the  next.  The  performance  indices  found 
in  Equations  (19)  and  (31)  will  not  change  with  incremental  control  inputs  -  the 
derivation  could  have  just  as  easily  began  by  replacing  u  with  Au  in  Equations  (9) 
and  (22). 

With  the  introduction  of  control  increments,  the  constraint  equation.  Equation 
(33),  must  be  modified  to  reflect  that  the  predictive  controller  will  be  determining 
a  control  increment,  not  an  absolute  control  input.  This  is  easily  accomplished  by 
substituting  the  equation  for  control  increments 


u{k  +  £)  =  Au{k  -I-  ^)  +  •  •  •  -I-  Au{k  +  q)-\-  u{k  -  1)  (36) 


into  Equation  (33).  After  simple  algebraic  manipulation  we  see  that  the  new  con¬ 
straint  equation,  involving  Au,  and  of  the  form  AU  <  6  is 


I 

-I 


7  0  0 

:  0 
I  ...  / 


Au{k) 


Au{k  -1-  q) 


^max 


^min 


I 

-I 


u{k  —  1) 
u(k  —  1) 


(37) 


3.5  Solution  Algorithm 

MATLAB  functions  were  developed  (see  Appendix  A)  to  solve  the  quadratic 
programs  given,  forming  the  predictive  controller.  SIMULINK  was  then  used  to 
couple  this  controller  with  the  system  model  to  complete  the  system.  Figure  1 
shows  the  complete  system  that  was  developed  to  solve  the  inverse  control  problem. 
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The  MPC-determined  input  is  passed  to  the  aircraft  model  where  the  actual  system 
output  is  determined.  This  output,  coupled  with  the  previous  input,  is  passed  to 
the  state  estimator,  and  the  process  begins  again  with  the  previous  inputs,  outputs, 
and  states  being  sent  to  the  controller  (MATLAB  fimction),  where  the  next  input 
to  perform  is  determined.  The  current  simulation  time  is  also  an  input  to  the  MPC 
so  the  controller  can  determine  where  in  time  along  the  setpoint  trajectory  the 
simulation  is,  and  the  setpoints  for  the  current  calculation  can  be  deduced  from 
that.  This  allows  for  the  prediction  horizon  used  in  the  MPC  calculations  to  be  both 
less  than  and  shifted  along  the  given  setpoint  trajectory. 

The  MATLAB  command  qp  was  used  to  solve  the  quadratic  problem  of  Equa¬ 
tions  (19)  and  (31)  once  the  state  prediction  matrices  were  calculated.  The  matrices 
constant  across  the  horizons  were  calculated  prior  to  executing  the  SIMULINK  pro¬ 
gram,  while  those  requiring  updates  across  time  were  calculated  within  the  MATLAB 
functions  (see  Appendix  A). 
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Once  the  MPC  SIMULINK  program  was  used  to  calculate  the  four  input  tra¬ 
jectories,  Se,ST,Sa,  and  Sr,  that  optimally  minimized  the  error  between  the  derived 
outputs  and  the  setpoint  trajectory,  a  second  SIMULINK  program  was  utilized  to 
calculate  all  of  the  outputs  which  were  contained  in  the  FDR  information  (see  Ap¬ 
pendix  B).  This  second  model  allows  inclusion  of  the  initial  conditions  of  the  flight 
since  the  model  used  in  the  MPC  program  was  initialized  to  steady  flight.  In  this 
manner,  the  inputs  required  to  drive  a  selection  of  the  outputs  to  the  setpoints  are 
determined,  while  the  effect  those  inputs  have  on  the  entire  set  of  known  output 
profiles  can  also  be  found. 
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IV.  Validation:  Distillation  Column 

This  section  presents  the  results  from  applying  the  inverse  model  predictive 
control  program  described  in  the  previous  sections  to  a  distillation  system.  This 
serves  as  a  validation  of  the  code  and  algorithm  developed  through  comparison  to 
results  foimd  under  a  previous  study. 

4.1  Scenario 

In  Morari’s  technical  notes  (9),  predictive  control,  based  on  a  step  response 
model,  is  performed  on  a  high-purity  distillation  column.  His  work  provides  a  unique 
opportunity  to  validate  the  results  obtained  from  the  algorithm  developed  for  this 
thesis.  The  algorithm  developed  for  this  thesis  will  be  applied  to  the  same  system, 
under  the  same  weighting  and  horizon  conditions,  to  validate  the  program  developed. 
Four  simulations  were  run,  each  varying  a  combination  of  the  input  weights,  the 
output  weights,  the  prediction  and  control  horizons,  and  the  setpoints. 

4.2  Model 

The  model  used  was  a  two-input,  two-output  transfer  function  system 

1 

75s  -f  1 

where  the  first  state  is  reflux  and  the  second  state  is  buildup.  This  thesis  utilizes 
state  space  models  so  the  transfer  function  model  given  was  transformed  into  a  state 
space  model  utilizing  MATLAB’s  tf2ss  command.  This  results  in  a  model  in  the 
form  of  Equation  (8),  for  continuous  time,  where 


-.0133  0 

1  0 

.0177  .0155 

- 1 

0 

0 

A  = 

B  = 

c  = 

D  = 

0  -.0133 

0  1 

.0144  .0146 

1 

0 

0 

1 _ 

0.878  0.864 
1.082  1.096 


(38) 
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4.3  Simulation 


The  first  set  of  results  determine  the  effect  of  the  control  weight  matrix,  R. 
These  effects  were  obtained  by  setting  the  setpoints  to  [1  0]^  throughout  the  horizon, 
9  =  5,p  =  20,  Q  =  /,  and  letting  R,  the  input  weight,  equal  0  and  I.  For  this  2-input, 
2-output  system  then. 


R  = 


0  0 
0  0 


,R  = 


1  0 
0  1 


(39) 


The  results  obtained  are  seen  in  Figure  2.  As  expected,  with  the  input  weight  set  to 
0,  the  controller  utilizes  large  control  inputs  to  drive  the  outputs  to  their  setpoints 
very  quickly;  while,  when  the  input  weight  is  set  to  I,  the  response  is  much  slower 
since  less  control  power  is  being  used. 


R-O  R-l 


Figure  2,  Distillation:  Effect  of  control  weight  matrix 


Next,  the  effects  of  the  output  weight  matrix,  Q,  are  determined.  Both  horizons 
and  the  setpoints  are  kept  the  same  as  the  previous  set  of  simulations,  R  is  set  to 
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the  identity,  and  two  values  are  used  to  find  the  output  weight  effects: 


Q  = 


100  0 
0  1 


,Q  = 


1  0 
0  100 


With  these  values,  the  validation  results  are  found  in  Figure  3.  From  this  figure,  it 
is  evident  that  the  output  associated  with  the  higher  weight  (100)  is  being  driven  to 
its  setpoint  much  faster  than  the  other  output,  which  has  a  weight  of  only  one. 


Figure  3.  Distillation:  Effect  of  output  weight  matrix 


With  the  weighting  matrices’  effects  found,  it  is  then  possible  to  determine  the 
effects  of  varying  the  setpoints.  This  is  done  using  four  different  setpoint  values,  r, 
remaining  constant  throughout  the  trajectory: 


1 

0 

0.88 

0.39 

r  = 

0 

,r  = 

1 

,r  = 

1.12 

,r  = 

0.59 

(41) 
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For  all  simulations,  q  =  5,p  =  20,Q  =  R  =  I.  These  results  are  found  in  Figures  4 
and  5.  Obviously,  as  the  setpoint  changes,  the  controller  drives  the  output  to  the 
desired  setpoint. 

From  these  runs,  it  is  evident  that  the  code  developed  is  valid.  The  results  are 
not  exactly  the  same  as  the  results  found  by  Morari,  since  the  exact  format  of  the 
performance  index  used  in  the  previous  study  is  not  known,  and  is  most  likely  not 
identical  to  the  index  used  in  this  thesis.  Therefore,  the  equation  being  minimized  is 
not  the  same,  and  slightly  differing  results  can  be  expected.  From  the  plots,  however, 
it  is  obvious  that  the  general  trends  obtained  from  the  algorithm  developed  are  the 
same  as  those  found  by  Morari, 

It  is  also  important  to  note  that  the  distillation  model  used  for  these  valida¬ 
tion  runs  is  very  different  from  the  typical  dynamic  model  for  an  aerospace  system. 
The  distillation  model  is  smaller  (two  states  versus  the  typical  eight  to  ten)  and  re¬ 
sponds  much  slower  than  expected  from  an  aircraft  model.  These  differences  aside, 
the  validation  runs  are  useful  for  two  reasons.  First,  as  stated  previously,  know¬ 
ing  the  results  of  a  completed  model  predictive  control  study  allowed  validation  of 
the  algorithm.  Second,  utilizing  this  model  shows  the  usefulness  of  applying  model 
predictive  control  to  a  variety  of  models. 
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Inputs  Outputs  ^  Inputs 


r-(1,0)  r-(0,1) 


4.  Distillation:  Effect  of  setpoint  changes  (r  =  [1,0]',  r  = 


r-(.88.1.12) 


Time  (sec) 


r-(.39,.59) 


Time  (sec) 


Figure  5.  Distillation:  Effect  of  setpoint  changes  (r  =  [.88, 1.12]', 


V.  Application:  Transport  Aircraft 

This  section  presents  the  results  from  applying  the  inverse  model  predictive 
control  program  to  a  transport  aircraft.  Two  actual  aircraft  mishaps,  Flight  427  and 
N827AX,  will  be  analyzed  to  find  the  degree  to  which  the  input  histories  can  be 
determined  from  the  FDR  data. 

5.1  Flight  427 

5.1.1  Scenario.  On  8  September,  1994,  a  Boeing  737-300  on  Flight  427 
crashed  while  on  approach  into  Pittsburgh  International  Airport.  The  aircraft  rolled, 
turned,  and  nosed  over  into  the  ground.  The  FDR  data  obtained  from  the  crash  con¬ 
tains  all  the  relevant  operating  conditions  of  the  aircraft  (altitude,  velocity,  heading, 
etc.)  but  none  of  the  control  surface  deflections. 

In  an  attempt  to  determine  these  inputs.  Parks,  Bach,  and  Shin  (12)  utilized 
the  Flight  427  information  and  an  assumed  set  of  control  inputs  to  reproduce  the 
outputs  that  were  provided  from  the  FDR  in  this  aircraft  accident.  Napolitano  (10) 
analyzed  the  same  aircraft  and  used  a  neural  network  simulator  and  “Virtual”  FDR 
to  reconstruct  control  surface  deflections. 

5.1.2  Model.  A  detailed  model  of  the  Boeing  737  is  not  available  in  open 
literature,  understandably,  due  to  its  proprietary  status.  However,  an  accurate  model 
estimation  can  be  made  from  aircraft  data  more  widely  available  in  the  current 
literature. 

The  FDR  data  obtained  for  this  particular  incident  suggested  that  the  aircraft 
was  in  a  steady  descent  at  about  6700  ft  100  seconds  before  impacting  the  ground. 
Therefore,  approximate  stability  derivatives  (13)  were  used  in  the  linearized  equa¬ 
tions  of  motion  to  provide  for  continuous-time  aircraft  dynamics,  where  in  state-space 
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form  (11) 


where 


X  =  Ax  +  Bu 


A  = 

Ajonj 

0 

B  = 

Blong 

0 

0 

^lat 

0 

Blat 

(42) 


(43) 


And  the  A  and  B  matrices  are  defined  in  terms  of  conventional  stability  derivatives 
(see  Appendix  C)  as 


Along  — 


-Ay  -Ay;  0  - Q 

Zn  Z'u)  XIq  0 

My  +  M^Zu  My;  +  M^Zy,  Mg  +  M^XIq  0 
0  0  10 


(44) 


^long  — 


A/ot  — 


Zse  Zsg. 

Mg^  +  MwZge  Mgg,  +  M.q,Zst 
0  0 

Yp/uo  Yp/uo  -{l-Yr/uo)  gcos{6o)luo  0 


Lp 

Np 

0 

0 


Np 

1 

0 


Blat  = 


Nr 

0 

1 

0  Ygjuo 
^Sa  Lg^ 
Nsy  Ng, 

0  0 

0  0 


0  0 
0  0 
0  0 
0  0 


(45) 


(46) 


(47) 
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The  aircraft  states  are 


where 


and  the  inputs  are 


=  [^Lg  xfatV 
Xiong  =  [A«  Aw  Aq  A0] 
xiat  —  [A/9  Ap  Ar  A(f>  AV»] 

Am  =  forward  velocity  change,  ftfs 
Aw  =  vertical  velocity  change,  ftfs 
Aq  =  pitch  rate  change,radls 
AB  =  pitch  angle  change,  rad 
Ap  =  sideslip  angle  change,  rad 
Ap  =  roll  rate  change,radl s 
Ar  =  yaw  rate  change,  radfs 
A<f>  =  roll  angle  change,  rad 
Arj)  =  yaw  angle  change,  rad 

U  =  [6e  St  Sa  SrY 


where 


Se  =  elevator  angle,  rad 
St  =  throttle  position,  percent 
Sa  =  aileron  angle,  rad 
Sr  =  rudder  angle,  rad 


(49) 


(50) 


(51) 


The  output  equation  y  =  Cx  can  then  be  formed  depending  on  which  states,  or 
linear  combinations  thereof,  axe  of  interest  later. 

Additionally,  the  following  equations  will  be  utilized  to  produce  outputs  in¬ 
cluded  in  the  FDR  that  are  not  simple  linear  combinations  of  the  states  described 
above.  The  absolute  velocity,  rate  of  climb/descent,  forward  and  vertical  accelera- 
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tions  are,  respectively, 


AVr  =  y/  Au^  +  At>2  +  Aw"^ 

Ah  =  UoA9  —  Aw 

(52) 

Ag^  -  Au/g 
Ag^  =  [Ah  +  i9)lg 

5.1.3  Simulation.  For  this  application,  both  rate  and  non- rate  performance 
indices  were  used  to  determine  the  input  histories  which  produced  the  data  contained 
in  the  FDR.  The  heading,  pitch,  and  roll  angles,  and  <j),  were  chosen  to  be 
included  in  the  performance  index  calculations  since  they  carry  most  of  the  attitude 
and  location  information  for  the  aircraft.  The  final  value  (at  t  =  100s)  was  extended 
for  just  enough  time  for  the  prediction  horizon  to  allow  for  one  final  input  calculation 
at  a  simulation  time  of  <  =  99s.  Before  being  included  in  the  setpoint  vector  passed 
to  the  MFC,  the  angles  were  resampled  to  the  same  discretization  interval  as  the 
plant  model,  AT =2  sec,  and  linearized  by  subtracting  the  initial  value  of  eaeh  angle 
throughout. 

5. 1,3.1  Non-rate  Performance  Index  Results.  Using  the  algorithm  as 
explained  in  Section  3.5,  and  with  the  following  parameter  values: 

L  =  observer  gain  matrix 

Q  =  diag{100  •  ■  •  100) 

R  =  diag{Q.l  •  •  •  0.1)  (53) 

p  =  20 

q  =  15 

the  four  inputs  (Se,ST,Sa,Sr)  that  optimally  drive  the  outputs  (^>,0,  ^)  to  the  set- 
points  are  found.  According  to  the  developed  algorithm,  these  inputs  are  next  used 
in  conjuction  with  the  second  SIMULINK  program  (see  Appendix  B)  in  order  to  de¬ 
termine  the  full  set  of  outputs  that  are  included  in  the  FDR  data.  These  calculated 
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outputs  are  shown,  along  with  their  associated  FDR  output,  in  Figures  6  through 

12. 

These  figures  show  the  good  agreement  between  the  outputs  obtained  from 
the  MPC-inverse  control  program  and  those  from  the  FDR.  The  average  deviation 
from  the  setpoints,  across  the  entire  trajectory,  is  only  0.4909,  0.0637,  and  0.1786 
radians  for  6,  and  (f>,  respectively.  However,  the  disadvantages  of  working  with  a 
linear  model  are  evident  in  these  plots.  Especially  when  analyzing  the  acceleration 
deviations  in  Figures  8  and  9  it  is  apparent  that  the  linear  model  used  here  is  having 
difficulty  matching  the  FDR  data.  In  fact.  Parks,  Bach,  and  Shin  (12)  found  that 
they  had  to  change  the  lift  and  drag  characteristics  of  the  aircraft  as  it  approached 
stall  in  order  to  match  this  data.  The  attempt  here  did  not  include  this  manipulation, 
without  significant  degradation  of  results.  The  poor  matching  with  the  accelerations 
is  due  to  the  lack  of  lift  and  drag  manipulation  in  this  study  and  due  to  uncertainty 
in  the  distance  from  the  measurement  accelerometer  (in  the  right  main  wheel  well) 
to  the  aircraft  center  of  gravity.  However,  overall,  the  match  between  the  results  and 
the  data  is  very  good. 

The  optimal  inputs  which  produced  the  above  output  results  are  shown  in 
Figures  13,  14,  15,  and  16.  Upon  analyzing  the  four  inputs,  elevator,  throttle, 
aileron,  and  rudder,  from  the  Flight  427  results,  it  appears  that  the  cause  of  deviation 
from  controlled  flight  was  a  rudder  hard-over  input.  This  is  consistent  with  several 
other  sources.  Most  notably,  Peirks,  Bach,  and  Shin  found  that  they  needed  to  input 
a  two  or  three-step  rudder  command,  much  as  seen  in  Figure  16,  in  order  to  obtain 
output  profiles  similar  to  the  FDR.  Moreover,  the  results  found  here  are  analytical 
-  not  assumed  input  trajectories  as  others  have  been. 
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Figure  7.  Output  Response:  Vy 
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Figure  8.  Output  Response: 
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Figure  9.  Output  Response: 
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Figure  10.  Output  Response:  ^ 


Figure  11.  Output  Response:  Q 
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Figure  12.  Output  Response:  <j) 


Figure  13.  Optimal  Input,  8^ 
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Figure  14.  Optimal  Input,  St 


Figure  15.  Optimal  Input,  Sa 
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Figure  16.  Optimal  Input,  6r 
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5. 1.3.2  Rate  Performance  Index  Results.  Using  the  same  solution 
algorithm  and  setpoints  as  in  the  previous  section,  but  replacing  the  performance 
index  with  the  terms  including  the  rate  term,  Ae,  in  the  MFC  the  results  of  the 
Flight  427  problem  can  be  determined  and  compared  to  the  results  from  the  MFC 
without  the  rate  term  included.  With  the  following  parameter  values  (Note  teririR 
distinguishes  parameters  used  in  the  rate-term  calculations): 

Lr  =  observer  gain  matrix 
Qr  =  diag{m  ■  •  •  100) 

Rr  =  diagitH.X  •  •  •  0.1) 

(54) 

E  =  d*O(gr(100  •  •  •  100) 

PR  =  15 
qR  =  15 

the  four  inputs  that  optimally  drive  the  outputs  (^,0,^)  to  the 

setpoints  are  found,  just  as  in  the  solution  for  the  non-rate  performance  index. 
Once  again,  these  inputs  are  used  in  conjuction  with  the  seven-output  program  to 
determine  the  full  aircraft  response.  These  calculated  outputs  are  shown,  along  with 
their  associated  FDR  output,  in  Figures  17  through  23.  Flots  of  the  four  optimal 
inputs  which  produced  these  outputs  are  in  Figures  24,  25,  26,  and  27. 

Only  small  differences  can  be  seen  between  these  plots  and  the  results  from 
the  previous  simulation.  However,  it  is  evident  that  there  is  improved  setpoint 
tracking  from  using  the  rate-included  performance  index.  Figures  21  through  23 
show  less  deviation  from  the  setpoint  than  Figures  10  through  12.  In  comparison 
to  the  non-rate  index  setpoint  deviations,  the  average  setpoint  deviations  across  the 
entire  trajectory  here  are  slightly  less.  The  following  table  details  these  differences. 
The  improvement  is  small,  however,  and  may  not  warrant  the  increased  calculation 
time  required  for  the  additional  terms.  (Run  times  on  a  Sun  platform  doubled  from 
approximately  45  to  90  seconds  with  inclusion  of  the  rate  term.) 
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Table  1.  Rate  eind  Non-rate  ] 

Deviation  Comparison 

Setpoint 

Rate  deviation 

Non-rate  deviation 

V’ 

0.4828 

0.4909 

e 

0.0626 

0.0637 

<!> 

0.1539 

0.1786 

The  lack  of  significantly  differing  results,  compared  to  the  non-rate  performance 
index,  may  be  due  to  the  fact  that  the  inclusion  of  the  rate  term  amounts  to  changing 
the  Q  weight.  The  introduction  of  the  error  weight  matrix,  E,  causes  a  direct  increase 
in  the  Q  weight,  and  its  inclusion  in  the  cross  terms  resulting  from  the  performance 
index  expansion,  along  with  the  time  shifts  between  the  T  and  as  well  as  the 
other  matrices,  introduce  what  is  equivalent  to  a  non-constant  Q  weight  across  time. 
So  slightly  better  results  (due  to  higher  weighting)  but  not  significantly  better  results 
can  be  expected. 


Figure  17.  Output  Response:  hn 
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Figure  20.  Output  Response: 


Figure  21.  Output  Response: 
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Figure  22.  Output  Response:  Or 
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Figure  23.  Output  Response:  <l>ji 
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Figure  24.  Optimal  Input,  Sej^ 
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Figure  25.  Optimal  Input,  ^Tn 
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Figure  26.  Optimal  Input, 


Figure  27.  Optimal  Input, 
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5.2  N827AX 


5.2.1  Scenario.  On  22  December,  1996,  an  Airborne  Express  Douglas  DC- 
8-63F  impacted  mountainous  terrain  at  about  3400  feet  in  Narrows,  Virginia.  The 
crew  was  performing  stall  maneuvers  on  a  test  flight  when,  after  one  of  the  planned 
stalls,  the  aircraft  plunged  unrecoverably  into  the  mountains.  As  with  Flight  427, 
the  FDR  data  obtained  from  the  crash  includes  all  of  the  operating  conditions  of 
the  flight  except  for  the  control  surface  deflections.  This  FDR  data  will  be  used  to 
determine  the  control  histories  during  the  flight  in  an  attempt  to  determine  what 
might  have  caused  the  incident. 

5.2.2  Model.  The  same  model  that  was  utilized  for  the  Flight  427  caJ- 
culations  will  be  incorporated  here.  Although  that  model  was  based  on  flight  at 
approximately  6700  feet  and  this  flight  is  at  approximately  13,500  feet,  the  flight 
characteristics  (and  therefore  stability  derivatives)  at  each  altitude  vary  only  slightly, 
so  minimal  differences  can  be  expected.  Additionally,  although  this  incident  involves 
a  Douglas  DC-8  and  the  previous  calculations  were  for  a  Boeing  737,  the  flight  char¬ 
acteristics  for  each  can  be  well  incorporated  into  a  general  transport  model  which  is 
being  utilized  here. 

5.2.3  Simulation.  Only  the  non-rate  performance  index  will  be  utilized  for 
this  flight  since  it  was  seen  from  previous  results  that  the  rate  performance  index 
did  not  provide  greatly  improved  results.  For  this  simulation,  however,  four  FDR 
outputs  will  be  used  as  the  setpoints:  the  three  flight  angles,  if),  9,  and  <j)  and  the 
altitude,  h.  As  before,  the  final  value  (at  t  =  119s)  was  extended  just  long  enough 
in  time  for  the  prediction  horizon  to  allow  for  a  final  calculation  at  a  simulation 
time  of  <  =  118s.  Also,  before  being  included  in  the  MFC  calculations,  the  angles 
were  resampled  to  match  the  discretization  interval  of  the  plant,  AT  =  2sec,  and 
linearized  by  subtracting  the  initial  value  of  each  setpoint  throughout  time.  It  should 
be  noted  that  while  the  downward  spike  in  the  altitude  plot  around  t  =  118s  is  not 
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feasible,  it  was  in  the  FDR  data  and  will  therefore  be  kept  in  the  data  for  the 
calculation. 

Using  the  algorithm  as  before,  and  with  the  following  parameter  values  (Note 
terniAx  distinguishes  terms  used  in  the  calculations  for  flight  N827AX): 

Lax  =  observer  gain  matrix 
Qax  =  diag{100  •  •  •  100) 

Rax  =  diag{0.1  •  •  •  0.1)  (55) 

Pax  =  15 
Qax  =  15 

the  four  inputs  {6e,ST,Sa,Sr)  that  optimally  drive  the  outputs  {tl),0,(f>,h)  to  the 
setpoints  are  found.  Using  these  optimal  inputs  as  the  inputs  to  the  seven-output 
SIMULINK  program  produces  the  outputs  plotted  in  Figures  28  through  34,  along 
with  the  associated  FDR  output. 

Again,  fairly  good  agreement  is  obtained  between  the  outputs  obtained  from 
the  MPC-inverse  control  program  and  the  FDR  data.  For  this  simulation,  average 
setpoint  deviations  across  the  entire  trajectory  are  0.2299,  0.3669,  and  0.2513  radians 
for  ip,  0,  and  <f>,  respectively,  and  868.9  feet  for  h.  While  the  agreement  is  not  as 
good  as  that  seen  from  the  Flight  427  results,  several  sources  of  error  are  entering 
the  problem.  As  stated  previously,  the  model  used  is  the  model  initially  developed 
for  a  Boeing  737  at  a  lower  altitude.  Stability  derivatives  more  representative  of  a 
DC-8  aircraft  at  altitude  would  produce  better  results.  These  differences  are  minor, 
but  are  contributing  factors  to  deviations  from  the  setpoints.  This  compounds  the 
linearization  factors  for  error  seen  in  the  Flight  427  results.  Here,  after  approximately 
60s,  the  aircraft  begins  to  develop  highly  oscillatory  flight.  At  these  conditions,  the 
linear  assumptions  at  the  root  of  the  model  begin  to  deteriorate.  It  is  at  this  stage  in 
the  simulation  that  the  outputs  depart  from  the  setpoints.  Plots  of  the  four  optimal 
inputs  which  produced  these  results  are  shown  in  Figures  35,  36,  37,  and  38. 
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Overall,  results  from  the  N827AX  simulation  show  good  setpoint  following. 
The  roll  data  is  matched  nearly  identically,  and  in  all  other  outputs,  the  trends  that 
exist  in  the  FDR  data  are  also  in  the  simulation  results.  Additionally,  Figure  35 
shows  a  stuck  elevator.  This  is  significant  since  a  stuck  elevator  is  the  suspected  cause 
of  the  incident  involving  N827AX.  The  results  seen  here  support  that  preliminary 
conclusion. 


Figure  28.  Output  Response,  Kax 
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Figure  30.  Output  Response,  Qxax 


Figure  34.  Output  Response,  <j>Ax 
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Figure  35.  Optimal  Input,  8ej^x 


Figure  36.  Optimal  Input, 
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Deter  (deg)  Detea  (deg) 


Figure  37.  Optimal  Input, 


Figure  38.  Optimal  Input, 
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VI.  Conclusions  and  Recommendations 

6.1  Conclusions 

This  thesis  proposed  that  a  Model  Predictive  Controller,  used  in  the  inverse 
sense,  might  prove  advantageous  for  use  in  aircraft  mishap  investigations.  As  a 
proof  of  concept,  the  FDR  data  obtained  from  Flights  427  and  N827AX,  in  which  the 
aircraft  crashed  for  unknown  reasons,  were  used  as  the  setpoints  for  the  MPC  and  the 
control  surface  inputs  that  may  have  caused  the  known  outputs  were  found.  Those 
inputs  resulted  in  flight  characteristics  showing  good  agreement  with  the  FDR  data. 
Average  deviations  from  the  setpoints,  across  the  entire  output  trajectory,  ranged 
from  0.06  to  0.5  radians  for  the  flight  angles.  Conclusions  then  are  two-fold:  first, 
MPC  appears  to  be  promising  for  future  analysis  of  aircraft  mishap  investigations, 
and  second,  analysis  of  these  flights  using  this  method  supports  previous  findings  for 
the  causes  of  the  incidents. 

Additionally,  when  the  results  from  the  MPC  using  non-rate  and  rate  terms 
included  in  the  performance  index  are  compared,  no  advantages  in  the  results  can 
be  seen  for  the  inclusion  of  a  rate  term  as  developed  for  this  application.  However, 
other  applications  may  warrant  the  use  of  a  rate  term  in  the  performance  index. 

6.2  Recommendations 

Two  logical  steps  for  exploration  into  this  area  of  using  MPC  in  the  inverse 
sense  become  apparent.  First,  it  might  prove  beneficial  to  replace  the  model  used  to 
find  the  actual  outputs  with  a  non-linear  model.  Some  deficencies  were  seen  when 
attempting  to  match  the  setpoints  to  the  FDR  data  which  might  be  eliminated  with 
the  aid  of  a  non-linear  model.  Second,  and  possibly  along  those  same  lines,  imple¬ 
menting  this  algorithm  with  a  high-performance  aircraft,  as  opposed  to  a  transport 
aircraft,  would  be  advantageous.  This  type  of  system  is  much  more  responsive,  more 
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dependent  upon  control  usage,  and  would  therefore  provide  an  additional  test  for 
MPC  in  the  inverse  sense. 

In  the  arena  of  including  a  rate  term  in  the  performance  index,  for  inverse 
control  or  forward  control  schemes,  additional  studies  need  to  be  performed.  In 
this  application,  its  inclusion  did  not  provide  a  dramatic  improvement  in  results. 
However,  before  a  definitive  statement  pro  or  con  on  the  rate  term  can  be  made,  it 
needs  to  be  used  with  other  applications.  In  the  case  of  a  more  responsive  aircraft,  it 
could  provide  improved  setpoint  tracking,  above  and  beyond  what  was  seen  in  this 
thesis. 
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Appendix  A .  MA  TLA  B  Functions 

Two  MATLAB  functions  were  developed  for  use  in  this  thesis.  The  first  is 
a  setup  file  which  calculates  all  time-constant  matrices  required  and  prepares  the 
setpoint  vector.  This  file  is  slightly  different  for  the  rate-inclusion  performance  index 
since  the  matrices  are  different  than  those  required  by  the  non-rate  performance 
index.  Both  are  listed  below: 


function  [G,H,F,MC,MQ,MS,P,Cona,Conb,p]=setupcon(A,B,C,L,Q,R,p,q,yd) 
y.  Setup. m  to  be  run  prior  to  running  the  Simulink  program 
y,  to  perform  the  inverse  problem.  Calculates  the  matrices 
'/  required  for  the  non-rate  performamce  index.  Of  the  form: 

y. 


y.  c  B 

y.  [(A-L*C)*B 
y.  G=[(A-L*C)-2*B 

y.  [  : 

y.  [(A-L*C)“(p-l)*B 

y. 


0  0  ] 

B  0  ] 

(A-L*C)*B  0  ]  (pn  X  qxi) 

:  0  ] 

(A-L*C)~(p-2)*B  ...  B] 


y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 


H=[L  (A-L*C)*L  (A-L*C)^2*L _ (A-L*C)“(p-l)*L]  '  (pn  x  xi) 

F=[(A-L*C)  (A-L*C)“2  ...  (A-L*C)“p]'  (pn  x  n) 

MC=I*C  (peta  x  pn) 

MQ=I*q  (peta  x  peta)  (Assumes  same  weights  across  time,  but) 

(q  can  have  varying  weights  on  certain  inputs) 

MR»I*R  (qxi  x  qxi)  (Assumes  same  weights  across  time,  but) 

(r  can  have  varying  weights  on  certain  inputs) 


y.  Constraints:  Cona(u)  <=  Conb 

%  Cona=I*Con  (p  x  qxi)  matrix  of  constraint  multipliers 
%  Conb=I*Max  (p  x  1)  matrix  of  constraint  maximums 

y. 


y. 

y. 

y. 

y. 


where  p  =  prediction  horizon,  q  =  control  horizon, 
n  =  #  states,  xi  =  #  inputs,  eta  =  #  outputs, 
q  =  output  weighting,  R  =  input  weighting. 


n=size(A) ;n=n(l) ; 
xissize(B) ;xi=xi(2) ; 
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eta=size(C) ;eta=eta(l) ; 


for  i=l :size(yd,l) 

ydt((i-l)*eta+l:i*eta,l)=yd(i, ; 

end 

yd=ydt ; 


ALC=A-L*C; 
for  i=l:p 

F((i-l)*n+l:i*n, :)=[ALC“i] ; 

H((i-l)*n+l:i*n, :)=CALC"(i-l)*L] ; 

MC((i-l)*eta+l :i*eta, (i-l)*n+l :i*n)=C; 

MQ((i-l)*eta+l:i*eta, (i-l)*eta+l: i*eta)=q; 
for  j=l:q 
if  i>=j 

G((i-l)*n+l:i*n, (j-l)*xi+l: j*xi)=CALC“(i-j)*B] ; 

end 

MR((j-l)*xi+l: j*xi, (j-l)*xi+l : j*xi)=R; 
end 

end 

P=2* (G ' *MC » *MQ*MC*G+MR) ; 

for  i=l:q 
for  j=l:q 
if  i>=j 

a((i-l)*xi+l:i*xi, (j-l)*xi+l: j*xi)=eye(xi,xi) ; 

end 

end 

end 

Cona= [eye (q*xi ,q*xi) ; -eye (q*xi , q*xi) ] *a ; 
save  datafile  P  F  H  MC  MQ  G  Cona  p  q  A  B  C  yd 

function  [G,H,F,MC,MQ,MS,P,Cona,Conb,p]=setupconrate(A,B,C,L,Q,R,E,p,q,yd) 
y.  Setup. m  to  be  run  prior  to  running  the  Simulink  program 
y,  to  perform  the  inverse  problem.  Calculates  the  matrices 
%  required  for  the  rate  performance  index.  Of  the  form: 

7. 

7.  C  B  0  0  ] 

7.  [(A-L*C)*B  B  0  ] 
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•/. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 

y. 


G=[(A-L*C)"2*B 
[  : 


0  ]  (pn  X  qxi) 

0  } 


(A-L*C) *6 

[(A-L*C)*(p-l)*B  (A-L*C)*(p-2)*B  . . .  B] 

H=[L  (A-L*C)*L  (A-L*C)“2*L _ (A-L*C) -(p-l)*L]  '  (pn  x  xi) 

F=[(A-L*C)  (A-L*C)“2  ...  (A-L*C)“p]»  (pn  x  n) 

MC=I*C  (peta  x  pn) 

MQ=I*Q  (peta  x  peta)  (Assumes  same  weights  across  time,  but) 

(q  can  have  varying  weights  on  certain  inputs) 

ME=I*E  (peta  x  peta) 

MR=I*R  (qxi  x  qxi)  (Assumes  same  weights  across  time,  but) 

(r  can  have  varying  weights  on  certain  inputs) 

Constraints:  Cona(u)  <=  Conb 

Cona=I*Con  (p  x  qxi)  matrix  of  constraint  multipliers 
Conb=I*Max  (p  x  1)  matrix  of  constraint  mciximums 


y. 

y. 

y. 

y. 


where  p  =  prediction  horizon,  q  =  control  horizon, 
n  =  #  states,  xi  =  #  inputs,  eta  =  #  outputs, 
Q  =  output  weighting,  R  =  input  weighting. 


n=size(A) ;n=n(l) ; 
xi=size(B) ;xi=xi(2) ; 
eta=size(C) ;eta=eta(l) ; 


for  i=l :size(yd,l) 

ydt((i-l)*eta+l:i*eta,l)=yd(i, :) ' ; 

end 

yd=ydt; 


ALC=A-L*C; 


for  i*l:p 

F((i-l)*n+l;i*n, :)=[ALC“i]  ; 
H((i-l)*n+l:i*n, :)=[ALC“(i-l)*L] ; 
MC((i-l)*eta+l:i*eta,(i-l)*n+l:i*n)=C; 
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MQ((i-l)*eta+l :i*eta, (i-l)*eta+l :i*eta)=Q; 

ME((i-l)*eta+l : i*eta, (i-l)*eta+l : i*eta)=E; 
for  j=l:q 
if  i>=j 

G((i-l)*n+l:i*n, (j-l)*xi+l : j*xi)=[ALC"(i-j)*B]  ; 

end 

MR((j-l)*xi+l : j*xi, (j-l)*xi+l : j*xi)=R; 

end 

end 

Gtil= [zeros (n,q*xi) ;G(1: (p-l)*n, :)] ; 

Ghat=G-Gtil; 

Ftil= [zeros (n,n) ;F(1: (p-l)*n, :)] ; 

Fhat®F-Ftil; 

Htil= [zeros (n,xi) ;H(1: (p-l)*n, :)] ; 

Hhat=H-Htil; 

P=2*  (Ghat » *MC '  *Mq*  (MC*Ghat+2*ME*MC*G)  +G » KMC '  *ME '  *Mq*ME*MC*G+MR)  ; 

for  i=l:q 
for  j=l:q 
if  i>=j 

a((i-l)*xi+l:i*xi, (j-l)*xi+l: j*xi)=eye(xi,xi) ; 

end 

end 

end 

Cona=[eye(q*xi,q*xi) ;-eye(q*xi,q*xi)]*a; 

save  datafile  P  F  Fhat  H  Hhat  MC  Mq  G  Ghat  ME  Cona  p  q  A  B  C  yd 

Within  the  SIMULINK  inverse  control  program  is  a  MATLAB  function  which 
performs  the  inverse  control  calculation.  Again,  it  is  slightly  different  for  each  index, 
due  to  the  differing  matrices.  These  inverse  problem  functions  are  listed  below: 

function  0UT=ip737con(IN) 

%  Performs  the  inverse  problem  (non-rate)  and  outputs  the  current 
y,  J  value  and  the  next  input  to  perform,  driving  the 
'/•  actual  output  to  the  desired  output. 

•/. 

y.  MS=[yd(k) '  yd(k+l) '  yd(k+2) '  . . .  yd(k+p) '] 
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•/. 


y.  Load  P  F  H  MC  MQ  G  Cona  p  q  A  B  C  yd  from  'datafile. m’ . 
load  datafile 

n=size(A) ;n=n(l) ; 
xi=size(B) ;xi=xi(2) ; 
eta=size(C) ; eta=eta(l) ; 

y.  Extract  u(k-l),y,x,t  from  input  vector  {u(k-l)  ,y,x;t} : 

Uk.l=IN(l:xi,l); 
y=IN(l+xi:xi+eta,l) ; 
x=IN(l+xi+eta:xi+eta+n, 1) ; 
t=IN(xi+eta+n+l,l) ; 

y  Calculate  where  we  are  along  the  desired  output  vector: 
y  Round  to  ensure  integer  value,  add  one  timestep  since 
y  t=0  is  row  1,  X  by  eta  to  step  along  yd  vector  for  pmax 
pmin=r ound (t/2)*eta+l; 
pmax=pmin+p*eta-l ; 

y  Store  the  current  yd  for  (Yd-y)  calculation: 
curyd=yd(pmin:pmin+eta-l) ; 

y  Calculate  the  MS  matrix  for  f  calculation 
MS=yd(pmin:pmax) ' ; 

f =2* ( (F*x+H*y) ' *MC ' -MS) ♦MQ*MC*G ; f =f ' ; 

Umax=[19.33*pi/180  69.5  20*pi/180  26*pi/180]'; 
Umin=-[21.83*pi/180  30.5  20*pi/180  26*pi/180] ' ; 
for  i=2:q 

Umax((i-l)*xi+l:i*xi,l)=[inf  inf  inf  inf]'; 
Umin((i-l)*xi+l:i*xi,l)=-[inf  inf  inf  inf]'; 

Ok.l ((i-l)*xi+l : i*xi , l)=Uk_l (1 :4) ; 
end 

Conb= ( [Umax ; -Umin] - [eye (q*xi , q*xi ) ; -eye (q*xi , q*xi ) ]  ♦Uk_ 1 )  ; 

du=qp(P,f ,Cona,Conb) ; 

J= . 5*du ' ♦P+du+f ' *du ; 
du=du(l:xi); 
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'/,  Establish  the  output  vector:  {curyd;du; J}: 

OUT=[curyd;du;  J] ; 

function  0UT=ip737conr(IN) 

*/,  Performs  the  inverse  problem  (rate)  and  outputs  the  current 
*/,  J  value  and  the  next  input  to  perform,  driving  the 
y,  actual  output  to  the  desired  output. 

y. 

y.  Y=  [yd (k)  '  yd (k+1) '  yd  (k+2)  '  . . .  yd (k+p) '  ] 

y. 

y.  Load  P  F  Fhat  H  Hhat  MC  MQ  G  Ghat  ME  Cona  Conb  p  A  B  C  yd 
y,  from  ’datafile. m' . 
load  datafile 

n=size(A) ;n=n(l) ; 
xi=size(B) ;xi=xi(2) ; 
eta=size(C) ;eta=eta(l) ; 

y.  Extract  u(k-l),y,x,t  from  input  vector  {uCk-l)  ;y;x;t} : 

Uk_l=IN(l:xi,l); 

y=IN(l+xi:xi+eta,l) ; 

x=IN(l+xi+eta:xi+eta+n,l)  ; 

t=IN(xi+eta+n+l,l) ; 

y,  Calculate  where  we  are  along  the  desired  output  vector: 
y.  Round  to  ensure  integer  value,  add  one  timestep  since 
*/,  t=0  is  row  1,  X  by  eta  to  step  along  yd  vector  for  pmax 

pmin=r ound ( t /2 ) ♦et a+ 1 ; 
pmzLX=pmin+p*eta-l ; 

y,  Store  the  current  yd  for  (Yd-y)  calculation: 
curyd=yd(pmin:pmin+eta-l) ; 

y,  Calculate  the  Y  matrices  for  f  calculation 
Y=yd(pmin:pmax) ; 
if  pmin-eta<0 

Ytil= [zeros (eta, 1) ;yd(pmin:pmax-eta)]  ; 
else 

Ytil=yd(pmin-eta:pmax-eta) ; 
end 
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Yhat=Ytil-Y; 

f =2*  (  (MC*Fhat*x+MC*Hhat  *y+Yhat )  '  ♦MQ*  (MC*G+ME*MC*G)  + 
(MC*F*x+MC*H*y-Y)  '♦ME’*MQ*(MC*Ghat+ME*MC*G))  ;f=f » ; 

Uiuax=[19.33*pi/180  69.5  20*pi/180  26*pi/180]’; 

Umin=-C21 .83*pi/180  30.5  20*pi/180  26*pi/180] ’ ; 
for  i=2:q 

UmeLx((i-l)*xi+l:i*xi,l)*Cinf  inf  inf  inf]’; 
Umin((i-l)*xi+l:i*xi,l)=-[inf  inf  inf  inf]’; 
Uk_l((i-l)*xi+l:i*xi,l)=Uk_l(l :4) ; 

end 

Conb=([Umax;-Umin]-[eye(q*xi,q*xi)  ;-eye(q*xi,q*xi)]*Uk_l) ; 

du=qp(P,f ,Cona,Conb) ; 

J= . 5*du ’ ♦P*du+f ’ *du ; 
du=du(l : 1+xi-l) ; 

y,  Establish  the  output  vector:  {curyd;du; J}: 

0UT=[curyd;du; J]  ; 
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Appendix  B.  Seven- Output  Program 

A  second  SIMULINK  diagram  was  utilized  in  the  solution  algorithm  in  order  to 
determine  the  full  seven  outputs  that  were  included  in  the  FDR  data.  This  program 
is  found  in  Figure  39  below. 


Figure  39.  Seven-Output  SIMULINK  Program 
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Appendix  C.  Additional  Transport  Model  Information 

The  following  tables  show  the  stability  derivative  values  used  to  develop  the 
transport  aircraft  model  (13). 


Derivative 

Value  (s 

Derivative 

Value  (s 

-.0292 

Ms^ 

-.4430 

-.2260 

Lp 

-3.190 

.1400 

Np 

.4990 

-.6740 

Lso 

3.840 

Yr 

0 

Ns. 

.4010 

Lr 

.9800 

LSr 

.3350 

Nr 

-.2150 

NSr 

-.3270 

Lr. 

-1.390 

Ysr 

.0250  ft 

Nr, 

Yp 

0 

Yp 

-31.50  ft 

Derivative 

Value 

-.0016  ft/s 

.894e-®  ft/s 

-.0007  1/ft 

Xs. 

Zs. 

Zst 

MSrp 

.8160e-« 

Additional  values  required  for  calculation  of  the  A  and  B  matrices  are: 

Uo  =  316.26  ftfs 
do  =  .0723  rad 
t  =  ^  ft 
g  =  32.2  /t/s^ 
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where  Uo  and  6o  are  taken  from  the  FDR  data,  g  is  gravity,  and  I  is  the  distance 
from  the  aircraft  center  of  gravity  to  the  location  of  the  measurement  accelerometer 
(estimated  here  to  be  zero  or  approximately  zero).  The  resulting  ten-state  A  and  B 
matrices  are: 

0  0  -1.0000  0  316.2612  0  0  0  0  0 

0  -0.0292  0.1400  0  -32.2000  0  0  0  0  0 

0  -0.2260  -0.6740  316.2612  0  0  0  0  0  0 

0  0.0002  -0.0011  -0.7097  0  0  0  0  0  0 

0  0  0  1.0000  0  0  0  0  0  0 

A  — 

0  0  0  0  0  -0.0996  0  -1.0000  0.1015  0 

0  0  0  0  0  -3.1900  -1.3900  0.9800  0  0 

0  0  0  0  0  0.4990  -0.1130  -0.2150  0  0 

0  0  0  0  0  0  1.0000  0  0  0 

0  0  0  0  0  0  0  1.0000  0  0 

(56) 

0  0  0  0 

0.4500  0.0003  0  0 

-4.9500  0.0000  0  0 

-0.4394  0.0000  0  0 

0  0  0  0 

0  0  0  0.0001 

0  0  3.8400  0.3350 

0  0  0.4010  -0.3270 

0  0  0  0 

0  0  0  0 

where  the  states  are 

=  [Ah  Au  Aw  Aq  AO  Av  Ap  Ar  A<l>  Ar/)]  (58) 
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as  defined  in  the  thesis  text. 


This  model  results  in  poles  in  the  following  table,  along  with  two  poles  at  the 
origin: 


Pole(s) 

Location 

Phugoid 

-.0150±.117h' 

Short  Period 

-.6915  ±  .6047i 

Spiral 

-.0142 

Roll 

-1.5891 

Dutch  Roll 

-.0506  ±  .9386z 
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